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ABSTRACT 

The correction of the influence of phase corrugation in the pupil plane is a fundamental issue 
in achieving high dynamic range imaging. In this paper, we investigate an instrumental setup 
which consists in applying interferometric techniques on a single telescope, by filtering and 
dividing the pupil with an array of single-mode fibers. We developed a new algorithm, which 
makes use of the fact that we have a redundant interferometric array, to completely disentangle 
the astronomical object from the atmospheric perturbations (phase and scintillation). This self- 
calibrating algorithm can also be applied to any - diluted or not - redundant interferometric 
setup. On an 8 meter telescope observing at a wavelength of 630 nm, our simulations show 
that a single mode pupil remapping system could achieve, at a few resolution elements from 
the central star, a raw dynamic range up to 10 6 ; depending on the brightness of the source. The 
self calibration algorithm proved to be very efficient, allowing image reconstruction of faint 
sources (mag =15) even though the signal-to-noise ratio of individual spatial frequencies are 
of the order of 0.1. We finally note that the instrument could be more sensitive by combining 
this setup with an adaptive optics system. The dynamic range would however be limited by 
the noise of the small, high frequency, displacements of the deformable mirror. 

Key words: Atmospheric effects - Instrumentation: adaptive optics - Techniques: high an- 
gular resolution - Techniques: interferometric - Stars: imaging - (Stars:) planetary systems 



1 INTRODUCTION 

The image obtained though a telescope is a convolution between 
the brightness distribution of the astrophysical object and the point 
spread function (PSF). In the Fourier domain, it is the multiplica- 
tion of the Fourier transform of the object and the Optical Trans- 
fer Function (OTF). To restore a correct image of the source, one 
therefore needs to know precisely the OTF. In the presence of static 
aberrations only, deconvolution is possible since the OTF can be 
obtained by observing an unresolved object. But when the OTF 
is changing with time - for example, in the presence of atmo- 
spheric turbulence -, calibration requires averaging the perturba- 
tions, whose parameters var y with time. Thi s is one of the reasons 
why speckle interferometry dLabevridl 19701) . one of the most well 
known post-processing techniques, still has some difficulty to cre- 
ate high dynamic range maps. 

This mainly explains why real-time adaptive optics (AO) sys- 
tems are a fundamental feature of large telescopes. With such sys- 
tems, the OTF of the telescope is controlled by the deformable mir- 
ror to be the same as the one of an uncorrupted telescope. However, 
technological limits appear for i) larger telescopes (e.g. extremely 
large telescopes), ii) shorter wavelengths (e.g. visible), or iii) ex- 
tremely high dynamic range imaging (extreme adaptive optics). In 



these three cases it may be advantageous to contemplate a comple- 
mentary approach using post-detection techniques. The combina- 
tion of both could be the solution to reach major scientific results 
like extra-solar planetary system imaging. However, to do so, such 
techniques would require the knowledge of the time varying OTF. 

In IPerrin et"al ] d2006h . we proposed a passive solution (i.e., 
requiring no real-time modification of the optical path) by using 
a remapping of the pupil. Single-mode fibers provide us with the 
technology allowing such a massive modification of the geome- 
try of the pupil, while keeping zero optical path differences. In ad- 
dition, they also provide perfect spatial filtering. Data collection 
and analysis are then similar to those u tilized for aperture masking 
dHaniff et al.ll 19871 : iTuthill et al.ll2000l) . with the noticeable advan- 
tage of having the flux of the whole entrance pupil, and the pos- 
sibility to completely disentangle instrumental from astrophysical 
information. 

In Sect.[2]we explain why imaging through turbulence is an ill- 
posed problem. After a recall of the principle of the instrument, we 
show in Sect.[3]how what was before an ill-posed problem can be- 
come a well-posed one. This translate into an algorithm described 
in Sect l3.3l Finally, we show in the simulations of Sect H] that we 
can therefore reconstruct perfect images with a dynamic range only 
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limited by detector and photon noise. In Sect. \5\ we conclude by 
giving a brief summary of our results. 



2 THE ILL-POSED PROBLEM OF IMAGING THROUGH 
TURBULENCE 

The image formed in the focal plane of a telescope is the convolu- 
tion of the object brightness distribution O(x) with the point spread 
function PSF(x) of the instrument: 



I(x) = 0(x) * PSF(x). 



(1) 



In the Fourier domain, the convolution operation is transformed 
into a multiplication, while the Fourier transform of the point 
spread function is the optical transfer function (OTF): 



li(u) = V(u) x OTF(w). 



(2) 



We choose ji{u) as the Fourier transform of the image, and V(u) 
as the Fourier transform of the object brightness distribution. This 
is to be in line with interferometric conventions, where it is also 
called the visibility function. The fact that the image depends on 
two unknown functions, V(u) and OTF(w), is the problem under- 
lying any image reconstruction algorithm; without adding further 
information, we have no way to disentangle the object from the 
PSF. 

Following an interferometric approach, we discretize the OTF 
to reduce the problem to a system of observables and unknowns. 
The OTF results from the autocorrelation of the complex values 
of the complex amplitude transmission inside the pupil. Thus, the 
OTF can be discretized by considering the pupil as being made of a 
number of coherent patches where phase and amplitude variations 
are negligible. Each patch is defined by a position vector r t and a 
complex amplitude transmission: 



G(r;) = g;e i( 



(3) 



with a phase 0; (e.g. atmospheric piston), an amplitude g t (e.g. scin- 
tillation) and where i 2 = f - 1. Each pair of patches selects 
one specific spatial frequency described by the frequency vector 
Uk - (?i - fj)/A; where r t and r 7 - are the location vectors of the 
patches projected in a plane perpendicular to the line of sight and 
A is the wavelength. Hence, the OTF at frequency vector u k is ob- 
tained by the relation: 



VTF(u k ) = Yj G{rdG{ rj )\ 



(4) 



where < B k is the set of aperture pairs which sample the &-th spatial 
frequency u k \ 



UiJ) : (r i -r j )/A = u k \ 



(5) 



This shows that the optical transfer function can be obtained from 
the knowledge of the complex amplitude transmission inside the 
pupil. Using Eq. (|2]), we can deduce a direct relation between the 
pupil transmission, the Fourier transform of the image, and the 
Fourier transform of the brightness distribution of the object: 



f i(u k ) = V(u k ))_ j G(r i )G(r j y. 



Or to simplify the notation: 



(6) 



(7) 



def def 

where, and hereinafter, we define: \i k = jd(u k ), V k = V(u k ), 
OTF, d = OTF(w,) and G t = G(r t ). 

The image reconstruction problem is then reduced to finding 
the unknowns {V,, G;; V/} given the //,'s. The ill-posedness of 
this task can be exhibited thanks to a simple example. In Fig.[T] the 
complex amplitude transmission in the pupil is binned into six dif- 
ferent elements (the G z 's). The autocorrelation of these six patches 
creates an OTF defined by a real value OTF (at central spatial fre- 
quency) and 9 complex values (OTFi, . . . , OTF 9 ). These OTF val- 
ues multiplied by the visibility function of the astronomical object 
(V =l, Vi,...,Vg) yield the Fourier transform of the observed im- 
age as one real value /i and 9 complex values (p\, . . . ,// 9 ). Since, 
by definition, the real value V is equal to 1 and since the phase 
of one of the complex amplitude transmissions can be arbitrarily 
chosen, the image reconstruction involves the computation of 29 
unknowns (15 complex values: G , . . . , G 5 , Vi, . . . , V 9 , minus an 
arbitrary phase) given only 19 measurements (the real value fi 
and the 9 complex values . . . ,/i 9 ). Our example demonstrates 
that the image reconstruction when the PSF is unknown is an ill- 
posed problem termed as blind deconvolution dThiebaut & Conanl 
Without adding further information, disentangling astro- 
nomical from instrumental information is impossible. 

To avoid having to disentangle the time-dependent OTF, a tra- 
ditional solution is to average its fluctuations. Over multiple obser- 
vations, the long exposure OTF is: 

OTF,- / Y, GiG l)' W 

To calibrate the OTF, the astronomer can observe a point-like star 
(i.e. such that V k = 1,V&) and apply the same averaging process. 
However, when phase variations become larger than wavelength, 
the average of the complex OTF tends t oward 0, and d econvolu- 
tion is impossible with a finite S/N ratio dThiebautll2005h . In prac- 
tice, long exposure images have a A/r effective resolution, where 
r 20 cm in the visible is Fried' s parameter. Two solutions have 
been proposed to overcome this problem and achieve the diffraction 
limit at A/D where D is the pupil diameter. The first solution is to 
correct the wavefront in real time so as to keep the wavefront pertur- 
bations smaller than the wavelength. This is achieved with an adap- 
tive opt ics system. The second solution, so called speckle interfer- 
ometry jLabevrid[l970k is to take short exposures with respect to 
the time scale of the perturbation, and to average the squared mod- 
ulus of the Fourier transform of the image. This way, the transfer 
function for the modulus of the Fourier transform of the observed 
brightness distribution becomes: 



OTF, 



(9) 



and is attenuated for spatial frequencies higher than r 1 A but dif- 
ferent from zero up to D/A. The Fourier phase of the observed 
brightness distribution can be ret rieved by mea ns of a third order 
technique such as the bispectrum dWeigeltll 1 977h . 

Here we propose an alternative approach. Instead of averaging 
the OTF, the goal is to have real time measurements of the complex 
amplitude transmission. Then, a post-detection algorithm can be 
used to obtain a calibrated OTF: 



OTF, 



^ GtG* 



(10) 



where G t and Gj are estimated complex amplitude transmissions. 
We have however demonstrated in this section that this information 
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Figure 1. This sketch illustrates Eqs. J4j and |7}. The OTF result from the autocorrelation of the pupil complex amplitude transmission, and the Fourier 
transform of the image is the multiplication of the OTF by the Fourier transform of the object observed. The unknowns are the 15 complex values 
[Go, • • • , Gs, V\, . . . , Vg}, whereas the observables provide only 9 complex values . . . ,jdg}. Deconvolution is therefore an ill-posed problem. 




Interferometric output 



Figure 2. In this instrument, the pupil (or an image of it) is subdivided into 
several sub-pupils whose outputs are injected into single-mode fibers. The 
fibers are then rearranged to create a new non-redundant pupil. Imaging on 
the detector is then performed as if no remapping had taken place. 



is unavailable on a simple image of the object. In ord er to recover 
the m issing information, we have proposed a system JPerrin et al.1 
l20Q6h . in which the telescope pupil is injected into an array of 
single-mode fibers, and rearranged into a new non-redundant exit 
pupil. 



3 FROM AN ILL-POSED TO A WELL-POSED PROBLEM 
3.1 The instrument 



The concept, proposed in IPerrin et is summarised in 

Fig. [2] Briefly, entrance sub-apertures collect independently the ra- 
diation from an astronomical source in the pupil of the telescope, 
and focus the light onto the input heads of single-mode optical 
fibers (of location vectors r t ). The radiation is then guided by the 
fibers down to a recombination unit, in which the beams are rear- 
ranged into a ID or 2D non-redundant configuration to form the 
exit pupil. Finally, the remapped output pupil is focused to form 
fringes in the focal plane where a different fringe pattern is obtained 
for every pair of sub-pupils. 



The amplitudes and phases of the fringes are measurements 
of the Fourier components given by the entrance baselines vectors 
(r t - rj)\ the Fourier components measured in the image are thus 
given by the relation: 



fitj = VkGiG*.. 



(11) 



where V k is the complex visibility of the observed object at the fre- 
quency u k - (ri-rj)/A, and G; and Gj are the complex transmission 
factors in the telescope pupil as defined in Eq. It is interesting 
to note the differences between this relation and Eq. J7}: thanks to 
the remapping, each measurement now corresponds to a single pair 
of sub-apertures. 

A second advantage of this setup comes from the fact that 
single-mode fibers act as spatial filters. As a consequence, the re- 
lation Gi = giQ 1(f>i is exact for each sub-pupil. Indeed, after being 
filtered by the fiber, the complex electric field, which is otherwise a 
continuous function, can be characterized by only two parameters: 
its phase and its amplitude. The discretization introduced in the pre- 
vious section is no longer an approximation which opens the way 
toward searching for an exact solution. 



3.2 On the unicity of the solution 

The fundamental idea of this paper comes from the fact that by 
using interferometric techniques, information can be retrieved to 
deconvolve an image from its PSF. As stated in Sect. |2] im- 
age restoration requires the knowledge of the complex transmis- 
sion terms Gj, which i s impossible with direct imaging. How- 
ever, iGreenawavl dl982h proved that the missing information can 
be en coded at higher frequencies (see also lArnotlll983l :l Arnot et al.1 
Il985h . Remapping enables an increase in the number of observ- 
ables ji while keeping the number of unknowns constant. This is 
possible since the complex visibilities V k only depend on the base- 
lines in the telescope entra nce pupil. They do not cha nge with a 
rearrangement of the pupil dTallon & Tallon-Boscll 1992b . 

This can be well understood in terms of unknowns and ob- 
servables. A remapped system is governed by Eq. dTTb . Providing 
M sub-apertures and R redundant entrance baselines, the number of 
complex unknowns are of M(M - l)/2 - R visibilities (V k terms), 
and M - 1 transmission factors (Gi terms). On the other hand, the 
number of measurements is M(M - 1)/2 (the iiq terms; with i ^ j). 
Hence, if R > M — 1, there are more observables than unknowns 
and the system of equations can be solved. 

The fact that both the V k and the G t terms can be deduced 
from the jiQ can be illustrated in a specific case. Fig. [3]is the same 
sketch as in Fig.Q] but with the 6 sub-pupils rearranged into a non- 
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Optical transfer function 



FT (Remapped object) 



FT(Image) 




Figure 3. This sketch illustrates Eqs. J4} and {TTJ in the case of a remapped pupil. As in Fig.[T] the OTF result from the autocorrelation of the pupil complex 
amplitude transmission, and the Fourier transform of the image is the multiplication of the OTF by the visibility values of the observed object. However, 
whereas there are still 15 unknown complex values {Go, . . . , G5, Vi, . . . , Vg}, the observables provide 15 complex values {//,... ,fi\5) and deconvolution is 
therefore possible. 



redundant configuration. This configuration was chosen to have the 
most compact configuration, but any other non- redundant co nfig- 
uration could have been used (see for example rGolavl[l97ll) . On 
the left panel are the complex transmission factors of the remapped 
pupil. The other panels show the Fourier transform values of, from 
left to right, the PSF, the astronomical object, and the image on the 
detector. The equation linking the observables fi it j to the unknowns 
Gi and V k is Eq. (fTTT ). Inversion of the resulting set of equations 
may be possible since the number of unknowns is larger than the 
number of measurements. This can be demonstrated by using the 
logarithm of terms in Eq. (fTTT ), which becomes 



InO/iyl) = ln(|V,|) + ln(g ; ) + ln(g,), 
for the real part, and 

= <txy t ) + - <P„ 



(12) 



(13) 



for the imaginary part. In these two equations, <I>() is the argument 
function, and g,, gj, 4>i and 4>j are as defined in Eq. ©. We obtain 
this way two sets of linear equations, one for for the phases: 



= M P • 

and one for the amplitudes: 

pn(N)] = M A ■ 



WV)] 



[hi(g)] 

mm)] 



(14) 



(15) 



where [ ] represents column vectors. M P and M A are two matrices 
containing 1, 0, and -1 values. Specifically to the example of Fig.|3] 
Eq. lff4l becomes: 
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The rank of this matrix is 12, while the number of unknowns is 15. 
The three terms of degeneracy are one for the absolute phase ref- 
erence, and two for the tip and tilt. Thus, by providing an arbitrary 



constrain on these three terms (the absolute phase is arbitrary and 
the tip and tilt only depend on the location of the image centroid), 
we can perform a singular value decomposition of the matrix and 
obtain from the measurements Hq a unique solution for the phase of 
the perturbations and object visibilities. The same method applies 
to the logarithm of the amplitude: 
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The rank of this matrix is 14, meaning all the amplitudes can be 
retrieved, except for the total brightness of the object. This param- 
eter can easily be constrained by normalizing the flux of the recon- 
structed image. The measurement of the amplitudes is an important 
issue since we have to correct for injection variability in the single- 
mode fibers. 

It is however clear that solving this system would require tak- 
ing the logarithm of the measurements. This would be very sensi- 
tive to additive noise. To get the best out of the data, it is better to 
fit the measurements using their complex values and Eq. (fTTT ). To 
do so, we developed a self-calibration algorithm which permits the 
use of thousands of snapshot all together to reconstruct an image 
up to the photon noise limit. 



3.3 A self-calibration algorithm for redundant arrays 

This section presents a self-calibration algorithm adapted to single- 
mode pupil remapping instruments, but also more generally to any 
kind of redundant interferometric array. Indeed, the equation /i/j = 
Vjfc Gi G* established in Sec 13.11 is common to all interferometric 
facilities. In the case of long baseline interferometry for example, 
fiij is the measurement of the complex coherence value between 
telescopes i and j, G t the complex transmission factor of telescope 
z, and the complex visibility of the astronomical object at the 
baseline formed by telescope i and j. 

The particularity of this self-calibration algorithm comes from 
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the fact that it gives complex visibility estimations without the need 
of a regularization term. This is possible thanks to the redundancy 
of the interferometric array. If one wants to make sure this algo- 
rithm is adapted to a specific interferometric facility, he would have 
first to establish the M P and M A matrices, and thus verify the unic- 
ity of the solution. 

In the next sections, we first start deriving an algorithm in the 
single-exposure case (Sect. 13.3.11 13.3.21 and I3.3.3K and then we 
show how to extend our algorithm to account for multiple expo- 
sures (Sect. 1533113331 . 



3. 3. 1 Log-likelihood 

Following the iGoodmanl Jl985h model for the noise of measured 
complex visibilities, we assume that different measured complex 
visibilities are uncorrelated and that, for a given measured complex 
visibility ///j, the real and imaginary parts are uncorrelated Gaus- 
sian random variables which have the same standard deviation. Un- 
der these assumptions and from Eq. (QjJ, the log-likelihood of the 
data is: 



transmissions: 



(16) 



k (i,f)€S k 



where G t and Gj are the complex transmissions for each sub- 
aperture and where & k is the set of sub-aperture pairs for which 
the interferences sample the k-th spatial frequency u k as defined by 
Eq. ([5]). In Eq. JT6K the statistical weights are: 



1 



1 



Var(ReOiy)) Var(Im( A / /J )) 



(17) 



Solving the image reconstruction problem, in the maximum 
likelihood sense, consists in seeking for the complex transmis- 
sions Gi and the object visibilities V k which minimize the value of 
€(V, G) given by Eq. fl5J. Unfortunately, the log-likelihood €(V, G) 
being a polynomial of 6th degree with respect to the unknowns (the 
Vfc's and the G z 's), proper means to minimize it have to be invented. 



3. 3. 2 Best object visibilities 

Given the complex transmissions G, the expression of £(V, G) in 
Eq. dT6t is quadratic with respect to the object complex visibilities 
V. Providing the complex transmissions G are known, obtaining 
the best object complex visibilities V is then a simple linear least- 
squares problem. The solution of this problem is found by solving: 



31 
dV k 



0, \fk 



(18) 



where, by linearity, the derivative of the real quantity £(V, G) with 
respect to the complex V k is defined as: 



d£ def d£ . d£ 
W k ~ dRe(V k ) +1 dlm(V k ) 



(19) 



Then: 

d£ 



— = 2 J] vv, 7 (g,G*^- R/ ) G*Gj 



= 2 V k J] W UJ \°i\ 2 \°j\ 2 ~ 2 Yj ^ G ^ G J ■ (20) 

Solving Eq. (TEl with the partial derivative expression in Eq. (f20t 
yields the best object visibilities given the data and the complex 



X w iJ G i G j»iJ 

X ^j\Gi\ 2 \Gj\ 2 



(21) 



Not surprisingly, this solution is a weighted sum of the complex 
visibilities measured by sub- aperture pairs which sample the &-th 
spatial frequency. 

The visibilities obtained by Eq. (1211 are not normalized. As- 
suming Vq corresponds to the null frequency, the following normal- 
ization steps insure that the sought visibilities are normalized: 

a = Vl (22) 

Vl <- Vila (23) 

Gj <- yfaGj. (24) 

It is worth noting that the likelihood remains the same after these 
re-normalization steps. 

3.3.3 Self-calibration stage 

Since, given the complex transmission factors, the best object com- 
plex visibilities can be uniquely derived, the initial optimization 
problem of £(V,G) can be reduced to a smaller problem which 
consists in finding the complex transmissions which minimize the 
partially optimized log-likelihood: 

^(G) d = *(V,G)| v=yt(G) (25) 

where V\G) is given by Eq. J2T}, possibly after the re- 
normalization steps. The second stage of our algorithm therefore 
consists in fitting the complex transmissions so as to minimize 
l\G) with respect to the complex transmissions. 

Since the criterion fl(G) is continuously differentiate, its par- 
tial derivatives cancel at any extremum of the criterion. Hence the 
so-called first order optimality condition that at the optimum of 
£\G) we must have: 

d£\G) 



= 0, V/. 



(26) 



Note that, unless £\G) is strictly convex with respect to the G/'s, 
Eq. i26t is a necessary condition but is not a sufficient one because 
it would be verified by all the extrema (local minima, local maxima 
or saddle points) of the criterion. 

Since depends on G, the chain rule must be applied to de- 
rive the partial derivative of l\G) with respect to the i-th complex 
transmission. For instance, the derivative with respect to the real 
part of the i-th complex transmission expands as: 

d^V, G) y d€(V,G) dRc(Vl) 

d Re(G/) v =v\G) ^ 



d£\G) 
8Re(Gi) 



dRe(V k ) 



F= yt (G ) dRe(G/) 



d€(Y, G) 



dlm(V k ) 



dlm(Vl) 
v=v\G) dRe(Gi) 



However, since minimizes £(V, G), we have: 
d€(V,G) 



dV k 



0, 



v=vHo 



and from the definition in Eq. (\9l of the partial derivative with 
respect to a complex variable, we deduce that: 



d€(Y,G) 



dRe(V k ) 



: and 



d€(V,G) 



dlm(V k ) 



0. 
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It follows that: 



d£\G) d€(V,G) 



d Re(G z ) d Re(G z ) F=F t (G) 

Since the same reasoning can be conducted for the derivative with 
respect to the imaginary part of the complex transmission and by 
definition of the derivation of a real quantity with respect to a com- 
plex variable given in Eq. JT9K the partial derivative of the partially 
optimized log-likelihood finally simplifies to: 

d£\G) d£(V,G) 



V=V\G) 



(27) 



In words, since V\G) minimizes £(V, G) with respect to V, the 
partial derivative of £\G) = £{V\G) with respect to G is simply 
the partial derivative of £(V, G) with respect to G into which the 
V is replaced (after derivation) by V\G). This property helps to 
simplify the calculations to come and, more importantly, shows that 
the global optimum must verify the modified first order optimality 
condition: 

di\G) d€(V,G) 



■ , V/ . 



(28) 



Finally, the partial derivative of 
transmissions G can be written: 



v=vHg) 

with respect to the complex 
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d€(V, G) 
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From this last expression, it is tempting to derive a simple iterative 
algorithm by solving Eq. ( |28l for G; assuming the other complex 
transmissions are known. The resulting recurrence equation 
is: 



G 



(n+l) . 



j:(i,j)e!B k 



f-UMSk 



(29) 

where Gj^ is the y-th complex transmission at n-ih iteration of the 

algorithm and is the k-th best object visibility computed by 
Eq. (f2~Tb with the complex transmissions estimated at n-th iteration: 



y(n) d £ f yt (Q(n)^ _ 



(30) 



3.3.4 Multiple exposure case 



Our algorithm can be generalized to the processing of multiple ex- 
posures of the same object. We assume that the instrument does not 
undergo any significant rotation with respect to the observed object 
so that the sampled spatial frequencies (the w^'s) and the corre- 
sponding sets of sub-aperture pairs (the £Vs) remain the same dur- 
ing the total observing time. We also assume that the object bright- 
ness distribution is stable so that the object complex visibilities (the 



V&'s) do not depend on the exposure time. At least because of the 
noise and of the turbulence, the measured complex visibilities and 
the instantaneous complex amplitude transmissions however do de- 



pend on me exposure ma ex t an a are respectively denoted fiij tt and 
Gi t . Under the boodmanl < fl985h approximation, the log-likelihood 
becomes: 

£(V, G) = YjTj Tj " Git G l Vk \ 2 (31) 

t k (iJ)eS k 

where the, possibly time dependent, statistical weights are: 
1 1 



(32) 

Var(ReO/^)) Var(Im(/i a ,)) " 

Since the object visibilities and the instrumental geometry do 
not depend on time, the same spatial frequency is measured at ev- 
ery exposure by a given pair of sub-apertures. Hence the condition 
given in Eq. (fl"8l) can be used to trivially obtain the best complex 
visibilities of the object: 



VI 



(33) 
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which simplifies to Eq. (I2TT) in the case of a single exposure. 

The updating formula for the time dependent complex trans- 
missions is obtained from the condition in Eq. J28l> by simply 
replacing the aperture index i by an aperture-time index i,t and 
straightforwardly : 



MUj)eS k 
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(34) 

which also simplifies to Eq. d29l) in the case of a single exposure. 
3.3.5 Algorithm summary 

Putting everything together, our algorithm consists in the following 
steps: 

(i) initialization: set n = and choose the starting complex 
transmissions G (0) ; 

(ii) compute the best object visibilities given the complex 
transmissions G (n) according to Eq. (1331) and, optionally, renormal- 
ize the unknowns; 

(iii) terminate if the algorithm converged; otherwise, proceed 
with next step; 

(iv) compute G ( " +1) by updating the complex transmissions ac- 
cording to Eq. d34l ); 

(v) let n := n + 1 and loop to step 2; 

Our iterative algorithm is very simple to implement and its 
modest memory requirements makes it possible to process over 
several thousands of snapshots all together. This is a requirement 
for faint objects or to achieve very high dynamic range. Yet, on a 
strict mathematical point of view, our algorithm may have a number 
of deficiencies. First, as already mentioned, the first order optimal- 
ity condition is necessary but not sufficient to insure that the global 
minimum (or even a local minimum) of l\G) has been reached. 
Other non-linear image reconstruction algorithms (blind deconvo- 
lution, optical interferometry imaging, ...) have the same restric- 
tion. In practice, checking that the algorithm converges toward a 
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similar solution for different initial conditions can be used to assert 
the effective robustness of the method. A second possible problem 
results from the updating of the complex transmission by Eq. ([341 . 
If a fixed point is discovered by the recursion, then it satisfies the 
necessary optimality condition; but it is also possible that the re- 
cursion gives rise to oscillations for the values of the sought pa- 
rameters. Note that our updating method is a non-linear one, anal- 
ogous to numerical methods for solvi ng linear equations such as 
the Jacobi and Gauss-Seidel methods (lBarrettetal.il 19941) . Unlike 
for our specific problem, it is however possible to prove the con- 
vergence of the recursion in the linear case. Again, the behavior 
of the algorithm in practice can effectively prove its ability to con- 
verge to a fixed point. If it appears that the update formula leads 
to oscillations, this problem can be completely solved by using an 
iterative optimization algorithm which guarantees that the partially 
optimized log-likelihood fl(G) is effectively reduced from one it- 
eration to another. Since the log- likelihood i s a sum of squares, a 
Levenberg-Marquar dt algorithm dMorell 19771) coupled with a trust 
region method dMore & SorenserJ 19831) would completely solve 
this problem. In practice, none of the numerous simulations we 
have conducted with our iterative algorithm have given rise to any 
of these convergence problems. 

Although derived in the specific case of a pupil remapping in- 
strument, our algorithm shares simi larities with the self-calibratio n 
method used in radio-astronomy dCornwell & Wilkinson! fl98lh . 
However, in our case, not only the phases of the complex trans- 
missions are miscalibrated and must be recovered but also the am- 
plitudes. Besides, we do not need a regularization term to over- 
come the sparsity of the (w, v) coverage by radio interferometers. 
Our derivation of the non-linear updating formul a in Eq. (f34t is also 
quite similar to the iterative method proposed bv lMatsonl ( [l99ll) for 
recovering t he Fourier phase from the bispectrum phase and later 
improved by Thiebaut (1994) to achieve better convergence capa- 
bilities. 



4 DYNAMIC RANGE ESTIMATIONS 

4.1 Analytical estimation of the photon noise limitations 

This system gives calibrated measurements of the spatial frequen- 
cies of the object. The advantage is straightforward. In classi- 
cal imaging, phase and amplitude errors create speckles in the 
image plane, therefore limiting the dynamic range. With a pupil 
remapped instrument, and assuming we are acquiring fast enough 
to freeze atmospheric turbulence, statistical errors due to photon 
and detector noise will t heoretically be the main limiting factor. 
iBaldwin & Haniffl ([2002 ) showed that the dynamic range of a re- 
constructed image is linked to the errors of the Fourier components: 



as for the phase (iGoodmanll 19851) : 



dyn 



(6V/V) 2 + (6(f>) 2 



(35) 



where n is the total number of data points, (SV/V) is the fractional 
error in amplitude, and Sep the phase error (in radians). For a total 
number of photons N ph and a number of apertures M, the amplitude 
of the fringe peaks in the Fourier transform of the image is equal 
to A^ h /M (assuming full coherence for the fringes). Considering a 
white photon noise of amplitude ^JN^, the signal-to-noise of the 
visibility modulus can be estimated with 



6V 
V ' 



(37) 



This leads to the following approximation of the dynamic range: 



dyn 



M(M - 1) 
2MWph 



(38) 
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Within these approximations, this result has the merit of clearly 
highlighting the advantage and the drawback of a single-mode 
remapping system: (i) an arbitrarily high dynamic range can be ob- 
tained anywhere in the image, providing the integration time is long 
enough; (ii) since additive noise is uniformly distributed on all the 
spatial frequencies, it is also evenly distributed across the whole 
field of view. To compare, an optical design isolating the photons 
of a bright object next to a faint companion - like a perfect adap- 
tive optics and coronographic system - would achieve a superior 
photon- wise dynamic range of dyn = Af ph . Extreme dynamic range 
imaging, as required for detecting extra- solar earths (dyn « 10 10 ), 
would therefore also require a long integration time with our sys- 
tem. 



4.2 Numerical simulations 

4. 2. 1 Simulation setup 

To perform these simulations, we used "YAO", an adaptive optics 
simulation software written by F. Rigaut using the Yorick language. 
This software allows us to generate corrugated wavefront with and 
without adaptive optics correction. In our simulations, the instru- 
mental setup corresponds to an 8 meter telescope under good see- 
ing condition (r « 20 cm at 630 nm) caused by four different lay- 
ers of turbulence at altitudes 0, 400, 6 000 and 9 000 meters. The 
wind speed ranges from 6 to 20 m/s depending on the altitude of 
the layer. The AO system is optimized to work in the near infrared. 
It consists in a classical Shack-Hartmann wavefront sensor and a 
12 x 12 actuator deformable mirror. The loop frequency has been 
set to 500 Hz, with a gain of 0.6 and a frame delay of 4 ms. The 
guide star is of magnitude 5. 

The remapping was done by dividing the 8 meter telescope 
pupil into 132 hexagonal sub-pupils of 80 centimeters in diameter 
each. They are filtered by the fundamental mode of single-mode 
fibers, coupled so as to maximize the injection throughput of an 
uncorrupted incoming wavefront. The injection efficiency in this 
case would be of 78%. However, at the operating wavelength of 
630 nm, the diameters of the sub-pupils are large compare to the 
Fried parameter (d/r « 4) and the coupling is expected to be much 
lower without adaptive optics (« 5% in these simulations). The 132 
sub-pupils are then rearranged in a non-redundant configuration, to 
produce a total of 8 646 sets of fringes on the detector. 

The total integration time was set to 40 seconds. However, 
because of the coherence time of the atmosphere, acquisition was 
done by sequences of short acquisition periods. We used a snapshot 
time of 4 milliseconds, during which we integrated the effects of 
phase variations. This was a way to account for fringe blurring due 
to dynamic piston effects. We also added to our measurements the 
photon noise as a Gaussian noise of variance the number of photons 
on each pixel. The number of photons was computed to account for 
a coupling efficiency of 5% into the fibers and a spectral bandpass 
of 60 nm. 

No chromatic fringe blurring was introduced since its influ- 
ence would highly depend on the chosen technical setup. More- 
over, there are several ways to avoid this problem. In the case of 
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Figure 4. Horizontal cuts of the point spread functions imaged in Fig. \5\ 
The source is a star of magnitude 5, observed at 630 nm. The acquisition 
setting consists of 10000 snapshots of 4 ms each on an 8 meter telescope 
(ro ~ 20 cm). The upper panel shows the point spread function after an 
uncorrected turbulence, while the lower panel shows the PSF after partial 
correction of the wavefront by an adaptive optics system. The three curves 
are obtained through the Eqs. {8}, J9j and fLOt . The last equation require 
a pupil remapping to get an estimation of the G; and Gj terms. At a few 
resolution elements from the star, AO systems and speckle techniques can- 
not achieve a dynamic range over 40. However, pupil remapping enables a 
perfect reconstruction of the PSF, with dynamic ranges over 10 3 . 

a ID non-redundant remapping, we recommend to spectrally dis- 
perse the fringes. If a 2D non-redundant reconfiguration is manda- 
tory due to a large number of sub-pupils, another solution could 
be to use a hyp er-chromatic magn ifier like a Wynne lens system 
(as proposed by iRibak et al.ll2004r) . At least, the use of a narrow 
spectral filter can completely avoid the chromatic blurring. 

4.2.2 Comparison with the speckle and adaptive optics 
techniques 

This first test was performed in order to demonstrate the recon- 
struction quality of the PSF. The astronomical object is a point-like 
source of Fourier transform values 1 (Vk = 1, Vk). The observ- 
ing wavelength is 630 nm. As described in the previous section, 
the simulated dataset consists of 10000 snapshots, each featuring 
8646 fringe sets. From each set of fringes, a complex coherence 
value fiij is extracted, and complex transmission factors G t are es- 
timated according to the iterative algorithm described in Sect. l3.3.3l 
Finally, we used Eq. dTOt to obtain the calibrated OTF, and thus the 
PSF. For comparison, the same corrugated wavefronts were used 
to obtain the PSF with two other techniques. The first one consists 
in averaging the complex instantaneous OTF by Eq. {§]). The sec- 



ond one (speckle interferometry) consists in averaging the squared 
modulus, removing photon noise bias and taking the square root 
as in Eq. ([9]). In this work, we did not make use of the bispectrum 
or closure phase since our object is point-like and therefore purely 
symmetric. The four left panels of Fig. \5\ represent the deduced 
PSF, with (lower panels) and without (upper panels) the use of an 
adaptive optics system. 

The first result confirm the usefulness of both speckle interfer- 
ometry and adaptive optics systems, even though we are observing 
at visible wavelengths. This is clearly seen in Fig. [4j where un- 
like uncorrected long exposure imaging, they permit the retrieval 
of spatial information at the diffraction limit of the telescope. Nev- 
ertheless, the Strehl ratio in images obtained by these techniques 
is very low and the background pollution remains important. With- 
out remapping, the best dynamic range achievable at visible wave- 
length is with a combination of AO and speckle interferometry 
technique, which provide a dynamic range of 40. This of course 
is in the case of a perfectly working near-infrared AO system with 
a 500 Hz frequency loop. Most current AO system are not used in 
that mode since any slight miss-calibration of the actuator influence 
function would make it useless at these wavelengths. 

On the other hand, a diffraction pattern is obtained by using 
our remapping instrument and our image reconstruction algorithm. 
The pattern differs slightly from a perfect Airy disc due to the 
hexagonal sampling of the OTF. At a few A/D from the central star, 
a dynamic range of 10 3 is obtained (see Fig. [4]). Using an adaptive 
optics system does not significantly modify the pattern, meaning 
that the dynamic range is limited by the shape of the perturbation 
free PSF. Our conclusion is therefore that the dynamic range could 
be further increased by choosing a different PSF. We did this in the 
following section over a complex astronomical object. 

4.2.3 Image reconstruction 

In this simulation, the astronomical object is a star surrounded by a 
protoplanetary disc. The disc has an exponential brightness distri- 
bution and a total flux of a hundredth of the star flux. Two compan- 
ions are also present, one with a flux of a thousandth, and the other 
of ten thousandths the flux of the star. 

As before, the G t complex transmissions are estimated by us- 
ing iteratively Eq. ( [29l . But unlike in Sect. 14.2.21 calibrated ob- 
ject visibilities are retrieved by our algorithm using Eq. (f33T ) from 
10000 short exposure images. We used these visibilities to recon- 
struct an image. But instead of doing a simple Fourier transform 
to obtain the equivalent of a dirty map, it is possible to choose 
the PSF according to a scientific goal. For example, an Airy disc 
PSF is good to achieve high angular resolution, while a Gaussian 
apodization is better for the dynamic range. Of course, any other 
visibility apodization can also be applied to obtain the most suit- 
able PSF. This system thus has the advantage of producing images 
clos e to what can be obtained with optic ally apodized pupil op- 
tics (iKuchner & Traurj|l2002l : lGuvonll2003b . while completely free 
of any atmospheric perturbations. This way, the dynamic range is 
no longer limited by the Airy rings of the diffraction pattern of the 
telescope. 

However, there is one perturbation this system cannot correct 
for; it is the photon noise. The phase corrugation being accounted 
for, the photon noise should be the theoretical limit of the dynamic 
range. To try to reach this limit, we used a Gaussian shape apodiza- 
tion on the visibilities. All of the image reconstructions in Fig. [6] 
are of the object described in the first paragraph of this section, 
but with different total brightnesses. The reconstructions from the 
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Figure 5. Point spread functions (PSF) of the instrument using different imaging techniques. The point-like object is a star of magnitude 5, observed at 630 nm. 
The acquisition setting consists of 10000 snapshots of 4 ms each on a 8 meter telescope. Left panels: average of the exposures over the total acquisition time 
(cf. Eq. (8j). Central panels: reconstruction obtained by using conventional speckle interferometry technique (cf. Eq. (9j). Right panels: PSF after remapping 
and correction by its estimation (cf. Eq. fLOt). Upper panels: the PSF is obtained trough an atmospheric turbulence of ro ~ 20 cm (D/ro ~ 40). Lower panels: 
PSF with the same corrugated wavefront but corrected by a simulated adaptive optics system. The field of view of each image is around 15 resolution elements 
(15A/D). The PSF quality goes from bad (long exposure without adaptive optics), to medium (AO correction and/or Speckle deconvolution), to perfect (after 
a remapping and a post-detection processing). 



Table 1. Dynamic range results 



Without AO 
Magnitude ^JN~J2 a D.R. b 



With AO 
V^W2 a D.R. b 






1.1 x 10 6 


0.9 x 10 6 


2.4 x 10 6 


1.8 x 10 4 


5 


1.1 x 10 5 


1.5 x 10 5 


2.4 x 10 5 


1.7 x 10 4 


10 


1.1 x 10 4 


1.3 xlO 4 


2.4 x 10 4 


1.6 x 10 4 


15 


1.1 x 10 3 


0.8 x 10 3 


2.4 x 10 3 


1.2 x 10 3 



a Theoretical dynamic range as predicted by Eq. d38V 
b Dynamic range obtained by taking the background rms of the recon- 
structed images of Fig. [6] 



left to the right are respectively of a central star of magnitude 15, 
10, 5, and 0. The upper set of panels are for a system without AO, 
while the four bottom panels are with AO activated. Fig. |7] gives 
a summary by plotting a horizontal cut of the different reconstruc- 
tions. Table Q] lists the dynamic ranges estimated on the images 
and compares them with the analytical approximation of the pho- 
ton noise established in Sect. 14.11 The dynamic range estimations 
are obtained by taking the inverse of the root mean square of the 
background of the image normalized by its pixel of maximum flux. 

When atmospheric turbulence is not corrected by an AO sys- 
tem, the reconstructions clearly highlights a dynamic range lim- 
ited by photon noise. For an object of magnitude 10, we achieved 
a dynamic range of 1.3 x 10 4 . This is comparable to what was 
calculated from the total photon count and the analytical relation 
dyn « This theoretical value, taking into account the 5% 

coupling efficiency in the fibers, was of 1.1 x 10 4 . A second inter- 
esting point is that when the brightness of the source increases by a 



factor 100 (delta mag = 5), the dynamic range increases by a factor 
around 10. This can be seen over 15 order of magnitude, indicating 
that the dynamic range has a roughly linear increase as a function 
of -y/Afph- This does 1) confirm the validity of the analytical estima- 
tion, and 2) prove the quality of our self calibration algorithm to 
restore the object visibilities. 

Noteworthy is the reconstruction of the system with a bright- 
ness of 15 mag. Due to a relatively low coupling efficiency in the 
fibers, an average of 250 photons are detected per 4 ms snapshot. 
According to Eq. d36l h it implies a S/N ratio of « 0.12 for each spa- 
tial frequency measurement. It is therefore impressive to note that 
even with such a vague knowledge of the the iterative algo- 
rithm of Sect. I3.3.3l is capable of restoring images with a dynamic 
range limited by the photon noise. This is possible since it com- 
bines the advantage of fitting the information from all snapshots 
together with a small, cpu friendly, recursive algorithm. 

When atmospheric turbulence is corrected by an AO system, 
the dynamic range does not increase linearly as a function of the 
brightness anymore. Indeed, in the lower panels of Fig.[6]and Fig.[7] 
the dynamic range is clearly limited by another factor. Thanks to the 
correction, the injection throughput is higher, around 23%. This al- 
lows an increase in dynamic range for the faintest source, but with a 
saturation around a few 10 4 . The reason for this unpredicted thresh- 
old appeared clearly when analyzing our simulations. During the 4 
ms integration time, mirror displacements happened twice to adapt 
for the atmospheric turbulence. These minor corrections , while in- 
creasing the spatial coherence (Cagigal&Canales 2000), have the 
drawback of introducing a phase noise overshoot at the frequency 
of the loop. This side effect results in a slight blurring of the fringes 
recorded on the detector and a bias on the measurements. To profit 
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Figure 7. Horizontal cuts of the reconstructed images showed in Fig. [6] 
The acquisition setting consists of 10000 snapshots of 4 ms each on an 8 
meter telescope (ro ~ 20 cm). The upper panel shows the reconstruction 
in the case of an uncorrected turbulence, while the lower panel shows the 
reconstruction with adaptive optics correction. On the upper panel, we can 
clearly see the dynamic range improving with the brightness of the source, 
from 10 3 (15 mag), to 10 6 (0 mag). At the highest dynamic range, we can 
see the exponential brightness decrease of the protoplanetary disc. Behind 
an adaptive optic, deconvolution is no longer limited by the photon noise, 
but by the high frequency differential piston introduced by the deformable 
mirror. 

from the advantages of the AO, solutions could either be to increase 
the acquisition rate, or, better, to adjust the control loop to decrease 
the amplitude of the overshoot. 



5 SUMMARY 

In this paper, we presented furth er investigation of the instrument 
introduced in lPerrin et all d2006h . 

• We established an analytical relation linking the Fourier com- 
ponents of the image, the Fourier components of the object, and 
the atmospheric perturbations (Eq. fTTh. We showed that inverting 
this equation allowed us to completely disentangle turbulence from 
astronomical information (Sect. 13. 2t . 

• We developed an analytical iterative self-calibration algorithm 
which enables inversion of the previously established equation over 
several thousands of snapshots simultaneously. This algorithm hap- 
pened to be robust, allowing visibility determination from Fourier 
component measurements with S/N ratio well below 1 (Sect. l3.3t . 

• Simulations of this system confirmed the validity of the algo- 
rithm and produced high dynamic range, diffraction limited, im- 
ages of complex astronomical objects. A dynamic range of the or- 



der of 10 6 was achieved at visible wavelength on an eight meter 
telescope and in the presence of good seeing conditions. Compared 
to actual AO systems, it represents an increase of around 10 4 in dy- 
namic range. We noted that the sensitivity of the instrument would 
increase by using an adaptive optics system, but at the price of a 
limitation in the achievable dynamic range (Sect.|4]>. 
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Figure 6. These simulations feature reconstructed image from the visibilities acquired from Eq. J2TV The acquisition settings consist of 10 000 snapshots of 
4 ms each on an 8 meter telescope in the visible (ro ~ 20 cm). The object is a central star with a protoplanetary disc and two companions, of relative flux a 
thousandth and ten thousandths (respectively, at the upper-right and upper-left of the central star - they are highlighted by the red circles on the upper right 
image). From left to right, the difference in the reconstructed images are due to the brightness of the object (magnitudes 15, 10, 5 and 0). The upper panels are 
reconstructed images in the case of an uncorrected wavefront (D/ro ~ 40). The lower panels are reconstructed images with an AO corrected wavefront. The 
field of view of each image is around 30 resolution elements (30 A/ D). The color scale on the right is linear, from to 3 x 10 -3 , and normalized to the flux of 
the central star. The reconstructions show dynamic range around or over 10 4 , except for the ones with a star of magnitude 15. In the three upper panels, we 
can clearly see the photon noise limit, evolving as a function of the brightness of the source. For faint sources, the lower-right panel show that dynamic range 
is increased by combining a remapping setup with an adaptive optics system. However, on a bright source, fringe smearing due to the high frequency control 
loop of the deformable mirror limit the dynamic range. An horizontal cut of these figures can be seen in Fig. [7] 
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